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Nuclei far from stability are studied by solving the Hartree-Fock-Bogoliubov (HFB) equations, 
which describe the self-consistent mean field theory with pairing interaction. Calculations for even- 
even nuclei are carried out on a two-dimensional axially symmetric lattice, in coordinate space. The 
quasiparticle continuum wavefunctions are considered for energies up to 60 MeV. Nuclei near the 
drip lines have a strong coupling between weakly bound states and the particle continuum. This 
method gives a proper description of the ground state properties of such nuclei. High accuracy is 
achieved by representing the operators and wavefunctions using the technique of basis-splines. The 
detailed representation of the HFB equations in cylindrical coordinates is discussed. Calculations 
of observables for nuclei near the neutron drip line are presented to demonstrate the reliability of 
the method. 
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I. INTRODUCTION 

The latest experimental developments as well as re- 
cent advances in computational physics have sparked re- 
newed interest in nuclear structure theory. In contrast 
to the well-understood behavior near the valley of stabil- 
ity, there are many open questions as we move towards 
the proton and neutron driplines and towards the lim- 
its in mass number (superheavy region). The neutron 
dripline represents mostly "terra incognita" . In these 
exotic regions of the nuclear chart, one expects to see 
several new phenomena 0,0 : near the neutron dripline, 
the neutron-matter distribution will be very diffuse and 
of large size giving rise to "neutron halos" and "neu- 
trons skins" . There are also expected collective modes 
associated with this neutron skin, e.g. the "scissors" vi- 
brational mode or the "pygmy" resonance. In proton-rich 
nuclei, we have recently seen both spherical and deformed 
proton emitters; this "proton radioactivity" is caused 
by the tunneling of weakly bound protons through the 
Coulomb barrier. With RIB facilities, nuclear theorists 
see an opportunity to study the effective N-N interaction 
at large isospin, as well as large pairing correlations. 

It is generally acknowledged that an accurate treat- 
ment of the pairing interaction is essential for describing 
exotic nuclei This work is specifically aimed to cal- 

culating ground state observables such as the total bind- 
ing energy, charge radii, proton and neutron densities, 
separation energies for neutrons and protons, and pair- 
ing gaps. There are several types of approaches in nuclear 
structure theory Q: for the lightest nuclei, ab-initio cal- 



culations (Green's function Monte Carlo, no-core shell 
model) based on the bare N-N interaction are possible 
yg. Medium- mass nuclei up to A ~ 60 may be treated 
in the large-scale shell model approach [6j. For heavier 
nuclei one utilizes either nonrelativistic |3, 0, or rela- 
tivistic 13, EH Ull mean field theories. The large pairing 
correlations near the driplines can no longer be described 
by a small residual interaction. It becomes necessary to 
treat the mean field and the pairing field in a single self- 
consistent theory. Furthermore, the outermost nucleons 
are weakly bound, which implies a large spatial extent, 
and they are strongly coupled to the particle continuum. 
These features represent major challenges for the mean 
field theories. We overcome these difficulties by solv- 
ing the HFB equations for deformed, axially symmetric 
even-even nuclei on a two-dimensional lattice, without 
any further approximations. So far, most of HFB cal- 
culations are based on spherical symmetry or up to a 
limited energy in the quasiparticle spectrum continuum. 
The importance of the axially symmetry lies on the abil- 
ity to emulate a big range of nuclei that are not spherical 
in nature, e.g. nuclei that have an non-trivial intrinsic 
deformation. We have developed and tested a new mean- 
field nuclear structure code that specifically addresses the 
computational challenges and opportunities presented by 
nuclei near the driplines. 

The present work represents an introduction of the 
splines method to the solution of the HFB approach in 
axial symmetry. For now, we will focus on the methodol- 
ogy of our approach. We outline here briefly the theoret- 
ical and computational details. We also present results 
for a few nuclear systems to demonstrate the convergence 
of the results. 



'Electronic address: edgar.teran@vanderbilt.edu 



^Electronic address: volker.e.oberacker@vanderbilt.edu 
t Electronic address: umar@compsci.cas.vanderbilt.edu 



II. STANDARD HFB FORMALISM 



2 



The many-body Hamiltonian in occupation number 
representation has the form 

h = J2 <i \ t U > e l £ j 

The general linear transformation from particle operators 
c,$ to quasiparticle operators (3,(3^ take the form [l2|: 

(#)-(£ ("> 

The HFB approximate ground state of the many-body 
system is defined as a vacuum with respect to quasipar- 
ticles 

h |*o > = . 

The basic building blocks of the theory are the density 
matrix 

Pij =< $o|c) c,|$ >= {V*V T ) ZJ , (2.3) 
and the pairing tensor 

iHj =< $o|c, c,:|$ >= (^*t/ T )y ■ (2.4) 
which give form to the generalized density matrix 

»-(-'*• i-V)- 

The HFB ground state energy including the constraint 
on the particle number N is given by 

E(K) =< $ \H- \N\$ > ■ (2.5) 

The equations of motion are derived from the variational 
principle 

5 [E{K) - tr A(K 2 - 11)] = , (2.6) 
which results in the standard HFB formulation 

[H,K] = , (2.7) 
with the generalized single-particle Hamiltonian 

«=("^> -<„%.)• (28) 

where h and A denote the mean field Hamiltonian and 
pairing potential, respectively, and the Lagrange multi- 
plier A is the Fermi energy of the system. 



A. Quasiparticle wavefunctions in coordinate space 

In practice, it is to convenient to transform the stan- 
dard HFB equations into a coordinate space representa- 
tion and solve the resulting differential equations on a 
lattice. For this purpose, we define two types of quasi- 
particle wavefunctions 4>\ and 4>2 , corresponding to each 
quasiparticle energy state E a : 

<fr{E a ,Ttrq) = J2U ia (2a) fair - aq), (2.9a) 

i 

ME a ,raq) = J2 V * a M™*) ■ ( 2 -9b) 

i 

The basis wavefunctions <\>i depend on the coordinate vec- 
tor r, the spin projection a = ±i and the isospin pro- 
jection q (q = +i corresponds to protons and q = — h to 
neutrons). 

The particle density matrix for the HFB ground state 
assumes a very simple mathematical structure in terms 
of 4>i and 4>2 H : 

p(raq,r'a'q') = < $ | ft(r'a'q') ijj(raq) |$ > 

= ^2 Pij M V(J <l) ^(r'cr'q') 

i,j 

oo 

= J2 ME a ,raq) <&(E a ,T> V 'q') . 

E a >0 

(2.10) 

The sum over the states E a has replaced the integral form 
of the equations, since the HFB continuous spectrum has 
been discretized for practical calculations. 
Instead of the standard antisymmetric pairing tensor k 
defined as 

K{raq,r'a'q') = < $ | ftr'a'q') tp{raq) |$ > (2.11) 

we introduce the pairing density matrix p which is Hermi- 
tian for a time-reversal invariant ground state and hence 
more convenient to use |J] : 

p{vaq,r' a' q') = {-2a') K{raq, r' - a'q') 

= (-2a') J2 ^3 M^q) 4>j{*' ~ *V) 

i,j 

oo 

= - ME a ,raq) (f>t(E a ,r'a'q') . 

E a >0 

(2.12) 

In principle, the sums go over all the energy states, but 
in practice a cutoff in the number of states is done up to 
a reasonable number (~ 60 MeV). 

Proceeding in analogy to the pairing density matrix, 
we replace the antisymmetric pairing potential A in Eq. 
(|2.8|) with the Hcrmitian pairing field h 

h{raq,r'a'q') = {-2a') A{raq, r' - a'q') . (2.13) 
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B. Normal density and pairing density 

From expressions (|2.10|l and 12.1211 for the density ma- 
trices we obtain the following expressions for the nor- 
mal density p q (r) and pairing density p q (r), which are 
defined as the spin-averaged diagonal elements of their 
correspondent matrices 

<7 

= EE^«M #,«(™*) . ( 2 - 14 ) 
a a 

Pq{ v ) = J2p(™<l,™<l) 

a 

= -EE^,«N) ^,aH • (2-15) 
a a 

The quasiparticle energy E a is denoted by index a for 
simplicity. The physical interpretation of p q has been 
discussed in the quantity [p 9 (r) AV/2] 2 gives the 
probability to find a correlated pair of nucleons with op- 
posite spin projection in the volume element AV . 

C. Kinetic and spin-orbit densities 

The kinetic energy density r q (r) is defined as a func- 
tional of wavefunctions (f>2 

r q {r) = V-VV 9 (r,r')|r=r' 



= EEl V ^, Q (rag)| 2 . 



(2.16) 



The spin-orbit density does not appear directly in the 
nuclear potential, but rather its divergence 

V • J,(r) = -i ^(v4 jQ (r, q)) ■ (V x a)4> 2 , a (r, q) . 

(2.17) 



D. Energy functional and mean fields 

Standard HFB theory yields the following expression 
for the total binding energy of the nucleus in its ground 
state, with contributions from the mean field and the 
pairing field 

EhFB =< §HFB\H\<&HFB >— E m f + E pair . 

To simplify the notation, we drop the isospin indices q, q' 
in this section and in the following section. In coordinate 
space, the mean field contribution is given by Q 

Emf = \\ d 3 r J dV^I^rVO + Mr^rV) ] 

<7,(T' 

x p(r'a',ra) , (2.18) 
and pairing energy contribution has the form 

E pair = \J d3r J d 3 r'Y,H^r'a') p(rV,ra) . 

(2-19) 

The quantity h denotes the mean field, i.e. the particle- 
hole (p-h) channel of the interaction 



fc(ro-,rV) = i(rcr,rV) + / d 3 r 2 J d 3 r' 2 ^ u (2) (rcr, r 2 cr 2 ; rV, r 2 '4) p{r 2 'a' 2 , r 2 a 2 ) 



(2.20) 



o-2,cr 2 

r 



(2) ~ 

where v l2 is the antisymmetrized two-body effective N- In a similar way, we find for the pairing mean field h, i.e. 

N interaction (see Appendix) . The kinetic energy matrix for the p-p and h-h channels of the interaction 
elements are given by 



t(ra,rV) = 5(r - r') [ -— V 2 

2m 



(2.21) 



h(ra, rV) 



d 3 r[ 



d 3 i 



E 



(2.22) 



(2) 

E. Pairing interaction assumes that the effective interaction v pair is local 

In practice, one tends to use different effective N-N w paU rcr > r ' ~ ^\ r i >(7 i> r 2 - = 5(ti-t) S a[ . a 
interactions for the p-h and for the p-p channel. If one x 5(r 2 ' — r') 5 IJ ^ a iVp{Ya,Y l — a') , 
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the pairing mean field Hamiltonian becomes 

fe(rcr,rV) = V p (r<r,r' -a') p(ra,r'a') . 

For the pairing interaction V p we utilize the form 

V p (ra,r' - a 1 ) = V 5{r - r') 5^ a , F(r) . 

This parameterization describes two primary pairing 
forces: a pure delta interaction (F = 1) that gives rise to 
volume pairing, and a density dependent delta interaction 
(DDDI) that gives rise to surface pairing. In the latter 
case, one uses the following phenomenological ansatz [2£j 
for the factor F 



and 



F(r) 



1 



Po 



(2.23) 



where p(r) is the mass density. 

The DDDI interaction generates the following pairing 
mean field for the two isospin orientations q = ±i 



1 



/ lg (m,rV) = - V w p q (r)F(r) J(r-r') S a , a , . (2.24) 

The pairing contribution to the nuclear binding energy 
is then 



E, 



'pair 



pazr ~ pair 



d A r 



T/(P) T/( n ) 



F(r) 



An important related quantity is the average pairing gap 
for protons and neutrons which is defined as [3J, |^ 

< A g > = -— trace (h q p q ^j 
= -jjr J d 3 r J dV^^rV) p q (r'a',ra) 

where N q denotes the number of protons or neutrons. In- 
serting the expression derived earlier for the mean pairing 
field we arrive at 



1 v o 

< A„ > = -- 



(?) 



2 N„ 



d 3 r p q (r) Pq (r) F(r) . (2.25) 



Note that the pairing gap is a positive quantity because 



V {q) < 0. 



F. HFB equations in coordinate space 

For certain types of effective interactions (e.g. Skyrme 
mean field and pairing delta-interactions) the particle 
Hamiltonian h and the pairing Hamiltonian h are diago- 
nal in isospin space and local in position space, 

h(raq,r'a'q') = 5 q>9 > S(r - r>« (r) (2.26a) 



h(raq,T'a'q') = 6 M , 6(r - r')h%,(r) . (2.26b) 

Inserting these into the above HFB equations results in 
a 4x4 structure in spin space: 



(h« - A) 
h q 



h" 



t>U \ _ F ( <t>l, a 



with 



(2.27) 



h q = 



hUr) h q L (r) 



h q = 



h q u (r) h q n (r) 



Because of the structural similarity between the Dirac 
equation and the HFB equation in coordinate space, we 
encounter here similar computational challenges: for ex- 
ample, the spectrum of quasiparticle energies E is un- 
bounded from above and below. The spectrum is dis- 
crete for \E\ < — A and continuous for \E\ > —A. For 
even-even nuclei it is customary to solve the HFB equa- 
tions with a positive quasiparticle energy spectrum +E a 
and consider all negative energy states as occupied in the 
HFB ground state. 



III. 2-D REDUCTION FOR AXIALLY 
SYMMETRIC SYSTEMS 

For simplicity, we assume that the HFB quasiparticle 
Hamiltonian is invariant under rotations R z around the 
z-axis, i.e. [H, R z ] = 0. Due to the axial symmetry of 
the problem, it is advantageous to introduce cylindrical 
coordinates {(f), r, z). It is possible to construct simulta- 
neous eigenfunctions of the generalized Hamiltonian 7i 
and the z-component of the angular momentum, j z 

T~L ^>n,n, q {^,r,z) = £„,n, g ip n ,n,q(<p,r,z) (3.1a) 
jz Tpn,a,, q {<i>,r,z) = Ml i/> n ,n,q(<i>>r,z) , (3.1b) 

with quantum numbers f2 = ±-|,±|,±|,... correspond- 
ing to each nth energy state. The simultaneous quasi- 
particle eigenfunctions take the form 



ip n ,n,q(4>,r,z) 




27T 



\e^f ^l q (r,z,i)J 



3.2) 



We introduce the following useful notation 



U 



(1,2) 
nflq 

'(1:2) 



(r,z) = <frnfi] q ( r > z >V , 
(r,z) = </>£$ q (r,z,l) . 



(3.3a) 
(3.3b) 
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From the vanishing commutator, [H,j z ], we can deter- 
mine the ^-dependence of the HFB quasiparticle Hamil- 
tonian and arrive at the following structure for the Hamil- 
tonian 



h(<f>,r, z) 



h' u (r,z) 



* h ' n (r,z) 



e + ^h' u {r,z) h' u (r,z) 
and the pairing Hamiltonian 



V e +% * h'tf (r, z) h' u (r, z) 



(3.4) 



(3.5) 



Inserting equations ({3.4(1 and 13.5J1 into the eigenvalue Eq. 
H2.27(l , we arrive at the reduced 2-D problem in cylindrical 
coordinates: 



/ (hi 



-ft 'v 


h 


Tl 












-A) 






h'll 




h> 


Ti -(/i' TT -A) 




-h' 

n U 


h'u 


h' 


U 


"it 








( U {1) \ 

1 U n.n,q } 






/ U n.n,q \ 




X 


r(l) 

n.Cl.q 
f/ (2) 

I r(2) j 


= E n> Q t g 




n,Q,q 

u {2) 

u n,U,q 
I r(2) J 



(3.6) 



Here, quantities /i', ft,', {/ and L are all functions of (r, z) 
only. Also, ft/ and ft' contain the implicit isospin depen- 
dence q. This is the main mathematical structure that 
we implement in computational calculations. For a given 
angular momentum projection quantum number O, we 
solve the eigenvalue problem to obtain energy eigenval- 
ues E nt Q tq and eigenvectors ipn,n,q for the corresponding 
HFB quasiparticle states. 



A. Representation of operators 

The Hartree-Fock Hamiltonian using the Skyrme ef- 
fective interaction can be written 



h = - V-- 

q 2m* 



V + U q + UcS q l 



2i 



(V • Iq + I g -V) - tB,- (V x a) . (3.7) 



where U q is the nuclear central field, Uc the Coulomb 
interaction, and the spin-orbit field part is given by 
B g - (V X a). The explicit form of these expressions for 
the case of the Skyrme interaction are included in the 
Appendix. Starting from the kinetic energy we apply the 
cylindrical form of the Laplacian operator to the stan- 
dard form of the wavefunction in Eq. I|3.2[l , to find 



hi o 

o t 2 2 



(3.8) 



whose elements are given by 



\ dr 2 r dr 



(0-1/2) 



dz 2 



+ 
f 



d£d_ d£_d_ 

dr dr dz dz 

\ dr 2 

dr dr dz dz 



(3.9a) 



ld_ _ ( (0 + 1/2) 
r dr \ r 



dz 2 



(3.9b) 



/ being the effective mass given in (|A.9ll . The local po- 
tential terms could also be cast into a matrix form 



where 



«ii 

«22 



V\\ = V 2 2 = U q + U C S q l 



(3.10) 



(3.11) 



Expressions for U q and Uc are given in the Appendix. 
The Hartree-Fock spin-orbit operator 



— lB q - (V x a) — > w q 
can similarly be written into the form 



Wn W12 
W21 W22 



(3.12) 



(3.13) 



with LLf 



Wl 2 



W 2 1 



O - 1/2 



r dr dz 

dr dz 



W22 



-B r 



r 

+ 1/2 



(3.14a) 
(3.14b) 

(3.14c) 

(3.14d) 



where B ri B z arc defined in the Appendix for the Skyrme 
force. 



B. Densities and currents 

Making use of the definitions for the normal density 
and pairing density, Eqs. (|2.14|) and l|2.15|l . we apply the 
bi-spinor structure of the quasiparticle wavefunctions to 
find the corresponding expressions in axial symmetry: 
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p q (r,z) = i- 2]T X £ [\U% q (r,z)\' Z + \L^ q (r.: 
V n>o ) E n >o 

P*(r,z) = ~ (2 jrf) xJ2 [U%(r,z) U%>* q (r,z)+L 
V n>o ) E n >o 



Lflqi 1 "' Z ) LnQ q (r, Z) 



(3.15) 
(3.16) 



Similarly, starting from definitions (|2.1ti[) and (|2.17() . we obtain expressions for the kinetic energy density and the 
divergence of the current 



n>o/ E n >o 



(SI - 1/2) 5 



dU. 



(2) 



dr 



U, 



(2) 



(fi + 1/2) 2 



dL 



(2) 

nQq 



(9r 



0*7, 



(2) 



0L 



(2) 



(2) 
2" 



(3.17) 



V-J g (r) = i(2 fi £Tg% 



\f)7"r( 2 ) art 2 ) 
uu nQ.q un 7in q 

dr dz 



dr dz 



nflq 



9. 



dU„ 



(2) 



dr 



017. 



1/2 (2) 



(2) 



nClq 



dz 

^ L nQq 

dr 



.(3.18) 



The total number of protons or neutrons is obtained by 
integrating their densities 



N q = 

with 

N n (l q = 



d 3 r Pq (v) = 2 £ £ N nQq 

\ Q>0 / £„>0 



(3.19) 



l^(»-,«)| a + l^(r,*)| a 



(3.20) 

which gives the contribution of the quasiparticle state 
\nflq > to the proton or neutron density. In the HF+BCS 



limit, N, 



J nQ.q- 



An analogous treatment of the 



pairing density yields 

Pa = 



\ = f d?r p q {v) =(2^) £ P ^ 
J \ Q>o ) E n >0 



(3.21) 



with the "pairing density spectral distribution" (with re- 
spect to energy and angular momentum) 



Pn^lq 



rdr 



dz 



U%(r,z)uif q (r,z) 



+ L%(r,z)L%(r )Z ) (3.22) 



In the HF+BCS limit, P nVtq -> (it • v) n n q . Finally, 
we state the normalization condition for the four-spinor 
quasiparticle wavefunctions as 



d 3 r ^ooW V-«o 9 (r) = 1 , 



(3.23) 



which leads to 



rdr 



dz 



\U^ q (r,z)\ 2 + \L%(r,z)\ 2 



+ \U%q(r,z)\ 2 + \L^ (r,z 



nVLq \ 



1 . (3.24) 



IV. LATTICE REPRESENTATION OF SPINOR 
WAVEFUNCTIONS AND HAMILTONIAN 

For axially symmetric nuclei, we diagonalize the HFB 
Hamiltonian (|3.tj|) separately for fixed isospin projection 
q and angular momentum quantum number SI. We solve 
the eigenvalue problem by direct diagonalization on a 
two-dimensional grid (r a ,zp), where a — l,...,N r and 
/3 = 1,...,N Z . (Because the grid usually extends from 
— z to +z, we have N r sw N z , so when referring to the 
number of points in the mesh we only mention the values 
of N r ). The four components of the spinor wavefunction 
ipir, z) are represented on the two-dimensional lattice by 
an expansion in Basis-Spline functions Bi(x) evaluated 
at the lattice support points. Further details about the 
Basis-Spline technique are given in Refs. [Il[ll. For 
the lattice representation of the Hamiltonian, we use a 
hybrid method 0, 0, 0] in which derivative operators 
are constructed using the Galerkin method; this amounts 
to a global error reduction. Local potentials are repre- 
sented by the Basis-Spline collocation method (local er- 
ror reduction). The lattice representation transforms the 
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differential operator equation into a matrix form 



N 



7C^" = <V« (n = l,...,N), (4.1) 



The HFB calculations are initialized using the density 
output from a prior HF+BCS run which results in fast 
convergence of the HFB code. Because the HFB problem 
is self-consistent we use an iterative method for the so- 
lution, and at every iteration the full HFB Hamiltonian 
is diagonalized. Typically 15-20 iterations are sufficient 
for convergence at the level of one part in 10 5 for the to- 
tal binding energy. The Fermi levels X q for protons and 
neutrons are calculated in every iteration by means of a 
simple root search using the equations Q 



/(A,) 
N n n q (\ q ) 



N q (\ q )-N q = 
2 E E N n n q (X q ) 



n>o/ e„>o 



(£n(lq — X q ) 



((£ n n g -X q r + Al Qq ) 



2 E n n q yjN nnq (l-N nnq ) , (4.2) 



where E denotes the quasiparticle energy, and £ is the 
equivalent single-particle energy (as defined by the BCS 
formalism). The quantity N in the last line of the equa- 
tion denotes the spectral norm of the density as defined 
in Eq. 13.201 The calculated value for X q is used in the 
next iteration. This process is repeated until convergence 
is achieved. 
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FIG. 1: Binding energy of O vs. energy cutoff. Top: cutoff 
in the quasiparticle spectrum, bottom: cutoff in the equiva- 
lent single particle spectrum (absolute value). All calculations 
were performed with B-Spline order M = 7, N r = 18 lattice 



points, angular momentum projection fl r , 
R = Wfm. 



and box size 



V. NUMERICAL PARAMETERS: 22 
CALCULATIONS 

In this section, we present a series of studies of the 
numerical parameters in axially symmetric HFB calcula- 
tions. In particular, we study the dependence of observ- 
ables on the equivalent single particle energy cut-off, the 
lattice box size, the number of mesh points, and the max- 
imum angular momentum quantum number £l m ax- The 
numerical tests are carried out for 22 0. This neutron- 
rich isotope has an N/Z ratio of 1.75 and is close to the 
experimentally confirmed dripline nucleus 24 0. 

A. Energy cutoff 

The numerical solution of the HFB equations on a 2- 
D lattice results in a set of quasiparticle wavefunctions 
and energies. The quasiparticle energy spectrum con- 
tains both bound and (discretized) continuum states. 
The number of eigenstates is determined by the di- 
mensionality of the discrete HFB Hamiltonian, which is 
N = (4 • N r ■ N z ) 2 , for fixed isospin projection q and an- 
gular momentum projection f2. In our calculations, we 



typically obtain quasiparticle energies up to about 1 GeV. 
It is well-known that zero-range pairing forces require a 
limited configuration space in the p — p channel because 
the interaction matrix elements decrease too slowly with 
excitation energy Q. One therefore introduces an en- 
ergy cut-off, either in the quasiparticle energy (E max ) or 
in the equivalent single particle energy (£ max ). Hence, 
in the case of zero-range pairing forces the infinite sum- 
mations over quasiparticle energies in the expressions for 
the densities /?, r, and current J are terminated at a 
maximum quasiparticle energy. 

The quantity E max has to be chosen such that the max- 
imum quasiparticle energy exceeds the depth of the mean 
field nuclear potential, and all of the bound states have to 
be included in the sums |3j . We follow the prescription of 
Refs. 0,113 to set the cutoff energy in terms of the equiv- 
alent single particle energy spectrum, £ n . For the Skyrme 
SLy4 force with pure delta-pairing, Dobaczewski et al. 
|31| deduced a pairing strength of Vq = —l70MeVfm 3 , 
with £ m ax = 60 MeV. We utilize the same parameters in 
all of our 2-D calculations. 

Even though £ max is a fixed parameter in the HFB 
calculations, it is interesting to analyze the sensitivity of 
observables to the value of the energy cutoff. In Fig. 2] 
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we plot the total nuclear binding energy for cutoff val- 
ues of Emax between 10 and 60 MeV and the same for 
E m ax from 20 to 60 MeV. We find that in both cases, 
the binding energy remains essentially constant for cutoff 
values of 40 MeV and above. Clearly, a cut-off below 40 
MeV results in significant changes in the binding energy 
because quasiparticle levels with large occupation prob- 
abilities are left out. This result is in agreement with the 
1-D radial calculations of Ref . 4] . 



increasing box size. Figure [5] shows that convergence is 
reached at R=10 fm. The behavior of the quasiparticle 
states with respect to the mesh boundaries has also been 
discussed in Ref. 0. For heavier systems, the box size 
has to be increased. A safe initial guess for R is about 
three times the classical nuclear radius. Tests also show 
that one may utilize the same mesh spacing for both light 
and heavy nuclei. 




6 8 10 12 14 



R(fm) 

FIG. 2: Bottom: Total binding energy of 22 as a function 
of the box size R. Top: Quasiparticle energies for states with 
large occupation probability (N n ) as a function of R. The 
spline order used was M = 9, N r = 19 grid points, £l max = f , 
and cutoff energy £„, = 60 MeV. 
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Number of Points (N ) 



FIG. 3: Total binding energy, Fermi level and pairing gap for 
neutrons in 22 vs. number of mesh points in radial direction, 
for fixed box size R — 8 fm. The quantities Q m ax and Emax 
are the same as in Fig. 2 



B. Lattice box size 

Using cylindrical coordinates, the lattice box size R 
defines the boundary in radial (r) direction; the box size 
in z direction is 2R. The value of R must be chosen large 
enough for the wavefunctions to vanish at the outer edges 
of the box and needs to be adjusted for optimal accuracy 
and computing time. Figure [5] shows the dependence of 
the binding energy on R for 22 0. The mesh spacing was 
kept at a constant value of Ar » 1/m. Figure [3 also 
presents some of the quasiparticle energy levels E^ with 
large occupation probability N n ; these levels correspond 
to low-lying states in the equivalent single-particle spec- 
trum. Evidently, the quasiparticle energies and the total 
binding energy converge in essentially the same way with 



C. Number of mesh points 

One of the major advantages of the B-Spline technique 
is that one can utilize a relatively coarse grid resulting 
in a lattice Hamiltonian matrix of low dimension. Fig. 
13 shows several observables as a function of the num- 
ber of radial mesh points, for a fixed box size R = 8 
fm. The binding energy, neutron Fermi level, and pair- 
ing gap for 22 O reach their asymptotic values at about 
18 grid points in radial direction. For the fixed (r, z) 
boundary conditions utilized in our work, the B-Spline 
lattice points show a (slightly) non-linear distribution, 
with more points in the vicinity of the boundaries. In 
the central region, the grid spacing for 18 radial points 
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is 0.75 fm. 



D. Projection of the angular momentum, 

It has been mentioned in the formalism section that 
all observables can be expressed by sums over positive 
j z quantum numbers fl > 0. The maximum value O mQX 
increases, in general, with the number of protons and 
neutrons (Z, N) and also depends on the nuclear defor- 
mation. There is no a priori criterion to fix £l max ; this 
numerical parameter needs to be determined from test 
calculations in various mass regions. We have performed 
calculations for 22 using Sl m ax values from 5/2 to 13/2. 
Figure ^displays the results for the total binding energy 
neutron Fermi energy and neutron pairing gap. All three 
observables converge at 11 = 9/2. 



a 2-D coordinate lattice by comparison with the 1-D 
coordinate space results of Dobaczewski et al. 0, 
for spherical nuclei. For this purpose we have chosen 
two very neutron-rich spherical nuclei: a light nucleus, 
22 0, with N/Z — 1.75 and a heavy system, 150 Sn, with 
N/Z = 2.0. Finally, we will also present results for 
a strongly deformed medium-heavy system, 102 Zr with 
N/Z = 1.55. This system was chosen because it allows us 
to compare our lattice results (which treat the continuum 
states accurately) to the "transformed harmonic oscilla- 
tor" (THO) expansion technique recently developed by 
Stoitsov et al. [23]. In this framework, a local-scaling 
point transformation of the spherical harmonic oscillator 
is used to expand the quasiparticle wavefunctions in a set 
of bound single-particle wavefunctions. 

A. Exotic spherical nuclei: 22 and 150 Sn 



1.5 




FIG. 4: Binding energy, neutron Fermi level, and average 
neutron pairing gap for 22 O vs. maximum angular momentum 
projection Q m ax ■ Box size R = 10 fm, N r = 18 and an energy 
cutoff of 60 MeV were used. 



VI. RESULTS 

In this section we present converged numerical results 
of our 2D-HFB code. Our main goal is to demonstrate 
the accuracy of our Basis-Spline expansion technique on 



In Table |U we compare our 2-D HFB results for the 
spherical isotope 22 O with the 1-D radial HFB method 
of Ref. 3]. Corresponding results in the 2-D THO ba- 
sis with 20 oscillator shells are also given. All calcu- 
lations were performed with the Skyrme SLy4 force in 
the p-h channel and a pure delta interaction (strength 
Vq = —170MeVfm 3 ) in the p-p channel, corresponding 
to volume pairing. The table lists several observables: 
the total binding energy (for comparison, the experimen- 
tal value is — 162.03MeV), the Fermi level for protons 
and neutrons, the neutron energy gap (for protons, the 
gap is exactly zero in all three calculations), the rms ra- 
dius, and the quadrupole deformation (note that both 
2-D calculations predict essentially zero deformation). 
Overall, the results of the axially symmetric code of the 
present work agree with the other two calculations in all 
the observables. The binding energy predicted by our 
2D-lattice code is very close (within 40 keV) to the 1- 
D lattice result, while the THO method result differs by 
80 keV. The difference in the Fermi level for protons is 
due to different conventions in choosing this energy for 
magic numbers. We choose the Fermi energy to be the 



TABLE I: Calculations for 22 for HFB+SLy4. The axially 
symmetric calculations (2D) of this work used a box size R = 
10/m with maximum f2 = | and an energy cutoff of 60 MeV. 
The spherical calculation of Ref. iTjj was made with R = 
25 fm and a j = 4f • All calculations were made with a cutoff 
at 60 MeV. 

1-D [31] 2-D (THO) [32] 2-D (this work) 



B. E. (MeV) 


-164.60 


-164.52 


-164.64 


A n (MeV) 


-5.26 


-5.27 


-5.27 


A p (MeV) 


-18.88 


-18.85 


-18.16 


A„ (MeV) 


1.42 


1.41 


1.40 


A p (MeV) 


0.00 


0.00 


0.00 


Rrms (fm) 


2.92 


2.92 


2.92 


ft 


* 


0.00002 


0.0008 


E patr (n) (MeV) 


-2.85 


-2.78 


-2.75 
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TABLE II: Comparison of calculations for spherical nucleus 
150 Sri with HFB+SLy4. The 1-D calculations were made by 
Ref. |31| . using a box size R = 30 and a linear spacing of 
points of 0.25 fm, with j ma x of 4r. Calculations by the axially 
symmetric HFB 2-D code were made using a box size R = 
20 fm with iV r = 23, maximum Q = 4^. In both calculations 
the pairing strength Vo was set to -170 MeV fm? , and the 
energy cutoff to 60 MeV . 



1-D 



2-D 



B. E (MeV) 
A„ (MeV) 
A p (MeV) 

A, 



(MeV) 
(MeV) 



Rrms (fm) 

E pair (n) (MeV) 



-1129.0 
-0.96 
-17.54 
1.02 
0.00 
5.12 
* 

-10.452 



-1129.6 
-0.94 

-17.69 
1.00 
0.02 
5.13 
0.01 

-10.057 



mid point of the energy of the last occupied level and the 
first unoccupied level. 

We now present results for the tin isotope 150 Sn, a 
heavy nucleus far away from the valley of /3-stability 
which is close to the two- neutron drip-line. Table |H] 
gives a comparison of our 2-D results (which predict a 
very small quadrupole deformation 02 = 0.005) with 
Dobaczewski's 1-D radial HFB calculations. The box size 
used in the axially symmetric calculations was 20 fm in r 
direction and 40 fm in the z axis, whereas the 1-D code 
had a 30 fm radial box. Also, the density of points has a 
different meaning in the radial code, since it uses a differ- 
ent grid than the one used in the B-Splines technique for 
our 2-D code. For these calculations the resulting mesh 
spacing in the 1-D code was 0.25 fm, whereas the maxi- 
mum mesh spacing in the 2-D one was 1.1 fm. In the 2- 
D calculations an approximately 3000 x 3000 matrix was 
diagonalized for each and isospin value, and for each 
major HFB iteration. The full calculation required about 
30 HFB iterations. Like in the oxygen isotope, the agree- 
ment is very good; a possible source of small discrepancies 
is the fact that our 2-D code yields /?2 = 0.005 whereas 
the 1-D code assumes an exactly spherical shape. Table 
ITT1 also contains another interesting piece of information 
on 150 Sn: the neutron Fermi level A„ is located less than 
1 MeV below the continuum which shows the proximity 
of this nucleus to the two-neutron driplinc. 



B. Deformed neutron-rich nucleus: 



l Zr 



Our main motivation for developing an axially sym- 
metric code is to perform highly accurate calculations 
for deformed nuclei, including the continuum states. The 
zirconium isotope 102 Zr is a heavy nucleus with strong 
prolate quadrupole deformation in its ground state. Its 
neutron to proton ratio of N/Z = 1.55 places it into the 
neutron-rich domain although it is likely far away from 
the neutron dripline (in the 1-D spherical HFB+SkP ap- 



proximation [2J| the last bound nucleus in the chain is 
predicted to be 136 Zr). We have chosen this isotope 
primarily because our results can be compared to the 
stretched harmonic oscillator expansion (THO) method 
mentioned above which does not involve any continuum 
states. 

In Table ITTT1 we present the results of our 2-D HFB cal- 
culations in coordinate space with the results obtained 
by the THO method. A comparison of the total bind- 
ing energy of the system in both methods shows a dif- 
ference of about 21 keV which can be considered small 
in comparison to the absolute value of the energy (the 
experimental binding energy value is —863.7 MeV). The 
pairing strength parameter, Vo, used in each calculation 
also makes a difference. Other observables (Fermi levels, 
rms-radius and deformation fi^) agree quite well, also. 
However, we find differences in the energy gap values 
(A„, A p ); these may be attributed to the different den- 
sity of states used in the two methods see Eq. (|2.25|) or to 
the extrapolation of the oscillator parameter in the THO 
approach to deformed systems. 



VII. CONCLUSIONS 

In this paper, we have solved for the first time the HFB 
continuum problem in coordinate space for deformed nu- 
clei in two spatial dimensions without any approxima- 
tions. The novel feature of our new HFB code is that 
it takes into account high-energy continuum states with 
an equivalent single-particle energy of 60 MeV or more. 
In the past, this has only been possible in 1-D calcula- 
tions for spherical nuclei [4|. Current 3-D HFB codes in 
coordinate space, e.g. Ref. p?j. utilize an expansion of 
the quasiparticle wavefunctions in a truncated HF-basis 
which is limited to continuum states up to about 5 MeV 
of excitation energy. 

The Vanderbilt HFB code has been specifically de- 
signed to study ground state properties of deformed ax- 
ially symmetric even-even nuclei near the neutron and 



TABLE III: Comparison of calculations HFB + SLt/4 for 
la2 Zr with two different methods in the axial symmetry. The 
configurational space calculations (THO) were made by Ref. 
[32] with 20 oscillator shells and pairing strength of -187.10 
MeVfm 3 . Calculations by the coordinate space HFB 2-D 
code were made using a box size R — 12/m with N r = 19, 
maximum Q 
of 60 MeV. 



-y, Vo -170 MeV fm and the energy cutoff 



2-D (THO) 



2-D (this work) 



B. E. (MeV) 
A„ (MeV) 
Ap (MeV) 



A„ 
Ap 

ft 



(MeV) 
(MeV) 
(fm) 



-859.40 
-5.42 
-12.10 
0.56 
0.62 
4.58 
0.429 



-859.61 
-5.46 
-12.08 
0.31 
0.34 
4.58 
0.431 
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proton drip lines. The large pairing correlations near the 
driplines and the strong coupling to the continuum rep- 
resent major challenges for the numerical solution. We 
have solved the HFB problem on a two-dimensional grid 
in cylindrical coordinates (r, z) using a Basis-Spline rep- 
resentation of wavefunctions and operators. B-Splincs 
are a generalization of the well-known finite element tech- 
nique. By using B-Splines of order M = 9 (corresponding 
to polynomials of up to 8-th order) we are able to repre- 
sent derivative operators very accurately on a relatively 
coarse grid with a lattice spacing of about 0.8 fm. While 
our current 2-D lattices are linear, a major advantage of 
the B-Spline technique is that it can be extended to non- 
linear lattices [TR llfil ] which will be particularly useful 
for an accurate and efficient calculation of neutron skins 
in heavy nuclei. 

In this work, we have used the Skyrme (SLy4) effec- 
tive N-N interaction in the p-h channel, and a pure delta 
interaction (corresponding to volume pairing) in the p-p 
channel. We present results for binding energies, defor- 
mations, normal densities and pairing densities, Fermi 
levels, and pairing gaps. 

We have investigated the numerical convergence of sev- 
eral observables as a function of lattice box size, grid 
spacing, angular momentum Sl ma x, and we have stud- 
ied the sensitivity of the observables to the continuum 
cut-off. These test calculations were carried out for the 
neutron-rich isotope 22 O with N/Z =1.75 which is close 
to the dripline nucleus 24 0. 

Our HFB-2D code predicts a spherical shape for the 
neutron-rich nuclei 22 and Sn. In this case, our 
calculations can be compared with the 1-D radial HFB 
results of Dobaczewski et al. Q, and indeed there is 
good agreement between the two. We also present re- 
sults for a strongly deformed system, 102 Zr, in which 
case we present a comparison with the stretched oscilla- 
tor expansion method of Stoitsov et al. [2^, • 

We have implemented our code on an IBM-SP mas- 
sively parallel supercomputer. Parallelization is possible 
for different angular momentum states f2 and isospins 
(p/n). We will also study alternative numerical tech- 
niques, in particular damping methods which we have 
utilized for solving the Dirac equation on a 3-D lattice 
0. 

In the near future, we plan to investigate several iso- 
tope chains, with particular concentration on deformed 
nuclei. We also plan to study a variety of Skyrme pa- 
rameterizations for the mean field, and both volume and 
surface pairing. As more data from existing RIB facilities 
become available, it is likely that it will become neces- 
sary to develop new effective N-N interactions to describe 
these exotic nuclei. Furthermore, our 2-D HFB ground- 
state wavefunctions can be used as input into coordinate- 
space based QRPA calculations 33] to investigate collec- 
tive excited states of nuclei near the driplines. 
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APPENDIX: SKYRME PARAMETERIZATION 

The (density dependent) two-body effective N-N inter- 
action is given by 

vi 2 ) = t {l + x P a ) <y(ri -r 2 ) 

+ i ti (l + x x P a ) {«s(n - r 2 )k) + k 2 8{n - r 2 )} 
+ i a (l + x 2 P) k'-Jfa-ra) k 
+ I t 3 (l + x 3 P a ) p-< 5(ri-ra) 

+ i W (<7i + <t 2 ) • {k x 6(n - r 2 )k} , 

P being the exchange operator, and k, k relative momen- 
tum operators. This form of the interaction with parame- 
ters Xo,Xi,X2,X3,to,ti,t2,ts,t4,, has being changed to an 
equivalent one with b±, b[, 6 2 , b' 2 , b 3 , 63, 64, b' A , parameters 
poj . This is done through the transformation 



( tl \ 1 


' 4/3 


8/3 


-2/3 


- 


-2/3 


-4/3 


4/3 




4 


-8/3 


2 


V hx 2 ) \ 


v -2 


4/3 


-4 



-4/3 



(bi 

b 2 
\b> 2 



and 



(A.l) 



to 
toxo 

h 

hx 3 
U 



-bo - »b' 



-bo + TT^n 



16, 



3 
26 4 



-b' 



^b> 



2b' 



(A.2) 



The last equation only holds for certain forces, as shown 
in Table \N\ For forces like SKI and SKO 64 and b' 4 get 
different values. 



1. Energy density 

Calculation of the energy expectation value for an ar- 
bitrary interaction involves carrying out an integration 
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TABLE IV: Skyrme force parameters. Values for new para meters bo, b' , bi, b[, 62, b' 2 , 63, b' 3 , b 4 , and b' 4 have been calculated 
using relations IA,ll and lA,2l from old parameterization . |Tfll| . [20j . Numbers have been rounded to three decimal places. 



Force 


bo 


K 


61 


K 


62 


b'2 


63 


b' 3 


bi 


b'i 


SkM* 


-2764.025 


-1560.55 


68.75 


68.125 


170.625 


68.437 


3898.75 


1949.375 


65.0 


65.0 


z a 


-3145.945 


-3316.251 


64.495 


58.315 


148.877 


61.405 


5577.823 


6707.621 


61.845 


61.845 


SkT6 


-2145.863 


-1600.426 


0.0 


0.0 


110.25 


0.0 


4005.312 


3204.25 


53.5 


53.5 


SLy4 


-3526.790 


-3320.210 


32.484 


-49.289 


185.325 


62.665 


5776.007 


6385.639 


61.5 


61.5 


Skll 


1000.310 


869.809 


32.354 


-49.803 


-432.059 


-1136.719 


580.693 


-2810.714 


62.13 


62.13 


SkI3 


-2034.628 


-1424.936 


32.301 


-127.914 


100.074 


-124.799 


3336.309 


3632.793 


94.254 


0.0 


SkM 


-2231.708 


-1679.676 


32.271 


-75.310 


-121.462 


-528.369 


3814.977 


3991.101 


183.097 


-180.351 


SkP 


-3359.948 


-2322.346 


44.642 


89.284 


190.343 


140.223 


5100.600 


3185.341 


50.0 


50.0 


SkO 


-1882.032 


-608.585 


22.537 


15.075 


-72.754 


-358.023 


2660.027 


237.585 


176.578 


-198.749 


SkO' 


-2068.449 


-987.770 


19.156 


8.312 


41.250 


-128.648 


3132.384 


1192.344 


143.895 


-82.889 



over six dimensions in coordinate space. One of the pri- 
mary advantages of an interaction that contains a delta 
function, like the Skyrme one, is that the evaluation of 
such integral becomes substantially simplified, and it's 
reduced to a three-dimensional evaluation 

E = ($\H\<f>) = j d 3 rH(r) . (A.3) 

The Hamiltonian density Tt(r) is composed of several 
terms 



H = Ho + Hls + He 



(A.4) 



The kinetic energy and some of the density dependent 
terms in the Skyrme interaction are included in 

q 

1 

+ h (pr - j 2 ) - b[ ]T (p q r q - j 2 ) 
q 

TP^P + b j £/>,V 2 „, ■ (A.5) 



The current densities ( j , j q ) appearing in this term are 
identically zero for time independent states. The finite 
range spin-orbit terms have the form 

-Hls = - 64 pV ■ J - b' 4 <°<?( V ' J «) • ( A - 6 ) 

Q 

The Coulomb term contains an integral over the proton 
density as well as the Slater exchange term, 

1/3 

l 4/3 . (A.7) 



r 2 (I ) t/vM 1 



2. Single Particle Hamiltonian 



Several effective quantities appear in this equation. The 
effective mass is defined by 



2m* 2m 



+ h P ~ b l Pq 



the effective current density 

I, = - 2 bi j + 2 b[ j, , 
and the effective spin density 

B q = b' x J q + 64 Vp + 64 Vp q 



(A.9) 



(A.10) 



(A.ll) 



As previously indicated, all of the terms in Ea. (|A.10|) 
vanish for bound states. Also, the first term in Ea. (|A.ll|) 
is usually ignored. 

The effective nuclear potential for the Skyrme force is 
given by 



U q = b p - b' Q p q + bi t- b[ T q 

*p a - lJ £p 2 q + 2p a p q 
q 

64 V • J - 64 V ■ J«j 



y (« + 2 )^ +1 -| 



+ 6' 2 V 2 p q - 62 V 2 P 

and the Coulomb field is 



7 ~ 2 / j3 / Pp( r ') 2 

U c = e d r --e 

r — r 



1/3 



lp P (r)} 



1/3 



(A.12) 



(A.13) 



B r and B z from equations (|3 - 1 3|) for the spin-orbit part 
representation of the potential operator are given by: 



The Hartree-Fock Hamiltonian using the Skyrme ef- 
fective interaction can be written as 

h 2 

h q= - ; V + U q + U c -6 q ,p 

2m* 

+ i (V • I q + VV) - iB q - (V X a) . (A.8) 



B r = Bq • e r 
B z = B q ■ e z 



Vribip + b^pq) (A.14a) 
Vz{biP + b' iP q) , (A. 14b) 



64 and 64 values are shown in Tablc lTVI for different forces. 
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